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In this article, we develop an intuitive and efficient, numerical technique to solve the quantum 
evolution equation of generic lattice-refined models in loop quantum cosmology. As an application of 
this method, we extensively study the solutions of the recently introduced, lattice-refined anisotropic 
model of the Schwarzschild interior geometry. Our calculations suggest that the results obtained 
from the approach are accurate, robust and are in complete agreement with the expectations from 
the von Neumann stability analysis of the model. 
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o 

O I- INTRODUCTION 

(N 

^ Loop quantum cosmology (LQC) |I] is a symmetry-reduced application of loop quantum gravity |2j, a theory that 

O leads to space-time that is discrete at Planck scales. An important result of LQC is that it is free of singularities. 

, /^ i This is seen as a general consequence of the quantum evolution equation, which is a "difference equation" for the 

~f. wave function and does not break down in the deep Planck regime, where the classical singularity is located [3, . 

,__! Although, LQC models generically appear to exhibit good behavior on the Planck scale, some of these models 

run into problems in the semiclassical regime [H [10]. For example, in the semiclassical sector of the isotropic model 

'"""' with positive cosmological constant A, the Wheeler-DeWitt equation (which would be a close approximation to the 

^H difference equation of LQC) has solutions that oscillate on scales of a^/X, where a is the scale-factor. This scale becomes 

I smaller with the expanding universe, ultimately becoming small enough that it falls under the discreteness scale of 

l^r\ the LQC difference equation. Thus, such an LQC model would necessarily deviate from the correct semiclassical 

I— I behavior [10]. Similar issues arise in the context of anisotropic LQC models, such as Bianchi I LRS [6] and the 

Schwarzschild interior geometry [71, where the evidence of poor semiclassical behavior emerges in the form of a 

K^ generic instability in the solutions of the quantum evolution difference equation ^9 . 

vQ These issues relating to the semiclassical behavior of certain LQC models can be resolved by introducing lattice- 

QQ refinement |llj . The main concept of this modification is that the underlying lattice over which the quantum wave 

f^ function has support, should be refined during the evolution of the model. In general this leads to a new type of 

04 the quantum evolution equation, one that no longer has a uniform step-size, thus presenting new challenges for its 

I analysis and its solution. It should be emphasized that this modification is not necessary in the Planck regime, 

■f—i but only in semiclassical regime. This is because, near the classical singularity there are only a few action steps 

f — ■ of the Hamiltonian that are relevant, and therefore lattice-refinement has little impact. Research on solving and 

^D analyzing lattice-refined difference equations is of current importance because it would play an important role in 

further advancing our understanding of semiclassical LQC and also help with the development of effective equations 

and dynamics. In fact, preliminary results from lattice-refined models of LQC suggest that it plays an important role 

in the development of an LQC based candidate for dark energy [4 and also for the origin of cosmic infiation [5]. 

Based on full LQG, there is a heuristic way of understanding why lattice-refinement is essential. Consider the 
context of a cosmological model and consider a state of quantum geometry corresponding to a fiducial volume cell. 
As the scale-factor increases the number of vertices of the state (graph) would have to increase correspondingly. To 
compute the Hamiltonian constraint in the full theory, one would then need to calculate the holonomy around the 
faces of an elementary cube surrounding each vertex. Now, as the number of such elementary cubes contained in the 
specified cell grow with the universe, the area of their faces would decrease. And for this reason, the edge length 
which is used for the computation of the holonomy would decrease. Thus, the dependence of this edge length on the 
scale- factor can be thought of as an effect of the full theory on simpler cosmological models. This is the main idea 
behind lattice-refined models. 

In this work, we will present an intuitive and efficient numerical method for solving the Hamiltonian constraint of 
lattice-refined models in LQC. It should be noted that solving the Hamiltonian constraint with uniform discreteness 
structure i.e. without refinement, is relatively straightforward and there are a variety of analytic and numerical 
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techniques available [HIS]- With lattice-refinement, the naive iterative approach towards computing a solution, simply 
fails. The reason is as follows: due to the non-uniform stepping, computing the future value of the wave function 
requires you to know past of values of the wave function, that are different from the ones computed in the previous 
iteration! We will explain this in detail with an example in the next section. Due to this difficulty it is simply not 
possible to use the large variety of techniques available in the numerical methods literature for solving such difference 
equations. However, in our current method we will overcome this issue through the intuitive approach of using a local 
interpolation formula for the solution, that will enable us to compute precisely those values of the wave function that 
are needed, to push the evolution forward. In other words, we will use the values of the wave function that were 
computed in the previous iteration to obtain a local analytic approximation formula of the solution. Then we will use 
this formula to compute the values of the wave function that will be needed for the next iteration. In this manner, 
we can obtain a complete iterative solution to the lattice-refined quantum evolution equation for any LQC model. 

This article is organized as follows. In the next section, through an explicit example, we will elaborate on the 
challenge associated with solving a ID lattice-refined difference equation using a simple recursive scheme. We will 
then introduce our local interpolation based methodology and solve the same equation, effectively and efficiently. We 
will also solve the same relation using a different approach by performing a change of variable that will make the 
relation uniformly spaced (this can always be done in the context of ID relations, but not in general) and then make 
an explicit comparison of the two solutions. In particular, we will demonstrate that they agree very well in the regime 
of interest. Next, we will proceed to solve a 2D lattice-refined LQC Hamiltonian constraint - one that cannot be 
mapped over to an equi-spaced form - using the same local interpolation based approach. We will demonstrate that 
the results in that case, agree with our general expectations of the solution, in particular, with regard to its stability 
properties. All the equations mentioned above, that will be solved in this work are relevant to the lattice-refined 
LQC model of the Schwarzschild interior geometry. We will end this article with a discussion of our results and some 
remarks on related work |12j . 

II. SOLVING LATTICE-REFINED MODELS 

We will start this section by very briefly reviewing the lattice-refined LQC model of the Schwarzschild interior 
geometry. We will follow the notation and treatment as presented in reference [10 . Consider a lattice with N vertices 
that is adapted to the symmetry of the Schwarzschild interior geometry. If there are N^ vertices along the direction 
of the triad labeled by r and iV^ vertices in the spherical orbits of the symmetry group whose triad is labeled by fj,, 
then N — NrN^^ . Therefore, the step-size along /i would be 5/N^ and that along t would be 5/Nr- Following the 
process detailed in reference [TU] we arrive at the following expression for the Hamiltonian constraint for this LQC 
model with lattice-refinement: 



2^^^- 



^( V k + 2^A^r ^1 + VFl) (;^t.+25N-\r+25N-^ " ^ t.-25N-\r+25N-') 



+{^\r + 5Nr'\- ^\T-5Nr^\) ((m + 25iV-i)V;^^4,^-i,, - 2(1 + 2^H^N-^)ii^P^,r 
+ ((A.-25A^-i)W4.Ar-J 

+25N-\^\T-25Nr'\+ ^\) [^^^25N-\r~25N---^^.+25N-\r-2SN--) 

= 0. 

We now have to make further assumptions on how exactly the lattice spacing is changing with changing /i and r. The 
simplest case is that the number of vertices in each direction is proportional to the geometrical area of a transverse 
surface. This yields Nr oc \/\t\ and N^ ex yj/il. We will start with a detailed analysis of this case in this section. An 
alternate choice for N^ and N^ , which was also considered in ,1Q\ is based on the intuition that the number of vertices 
in a given direction is proportional to the geometric extension of that direction. The resulting difference equation is 
more complex, but has improved stability properties |10j and therefore it is much preferable as a model for further 
study. We will analyze and solve that version, later in this section. 

A. One-dimensional lattice-refined relations 



In the simpler context of Nr ex \/\t\ and Nfj_ ex vl/i], the above master difference equation is variable separable, 
which makes it possible to rewrite it in the form of two \D difference equations (one with /^ as the independent 



variable and the other with r) that have non-uniform stepping. If we write the wave function ip^j,.r — Af^B-r then these 
two ordinary difference equations take the form shown below: 



(1) 



{^. + 25N-^)A^^,,^-^ - 2(1 + 2^H^N-^)^iA, + {^,- 25N-^)A^_,,^-^ = A(A^+,,^-. - A^^_,,^-^)25N-^ (2) 

Here A is the separation constant. Now, consider the Br relation that appears above in equation (llj). In order to 
solve it, we can try using a naive recursive method, but we will demonstrate here that such an approach will be 
unsuccessful. 
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FIG. 1: We plot the solution of the Br relation for different values of the separation constant A, using the interpolation method 
as presented m the text. The dashed-line corresponds to a value of X — 1.0, the dotted line corresponds to X = 2.1 and the 
solid-line is for X — —5.1. It is interesting to note that as the value of t increases, all the plots start to appear like a solid and 
continuous line. This is because the density of the points on each plot increases with increasing r (the stepping in r decreases 
as S/Nr). 



To keep the discussion simple, assume 2S = I. Lets start with some initial data for Br'. B 



i-i/\/T 



= Bq = and 



Bi = 1. We have enough data to find the third point that is present in the Br relation (HI i.e., B-^_^^,^ 
parameter t can advance only in steps of Xj \/t. Therefore, the next value that t can take is r-|- l/y'r 



With T = 2, the Br relation has terms like B. 



2-1/V2' 



Bi and B, 



i^-vl 



= B2. Now, the 
= 1 + 1/71=2. 
2- We just computed the value of B2, but to use 

However, there is 



the relation above, in order to evaluate B^ii,^ we would now need to know the value of -B2-1/ 
no direct way of determining this value from the already known data i.e., Bq, Bi and B2. Thus, the recursive method 
of solving this difference equation seems to fail because the previous iteration is unable to provide us with data points 
needed for the next iteration. We resolve this problem in our method by performing a local interpolation to obtain an 
analytic formula using the values we already know, including the one we just computed from the relation itself. Then, 
we evaluate this interpolation formula at the points where the data values are needed in order to progress the evolution 
forward. More specifically, we proceed as follows: A least-squares fit is done with the points Bq, Bi and B2 to obtain 
a locally accurate formula for Br. This is then used to evaluate B2_i/^ and thus, ^2 , j^ ; /j can be evaluated. Now, 
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FIG. 2: We plot the percent- difference d, in the solutions of the B^ relation for different types of interpolation schemes. The 
solid-line is the percent- difference between the solution based on the least-squares fit that we usually use, with the solution 
computed using a second-degree polynomial fit. The dotted-line is the percent- difference between the solution based on the least- 
squares fit, with the solution computed using a spline interpolation. The differences between these solutions are at the 0.1% 
scale. This is a clear indication of our interpolation based method being very robust. 



this computed value is used to get a revised fit formula, that is accurate in the local neighborhood of r = 2 — l/\/2, 
2 and 2 + l/\/2. Then, this revised formula is used to evaluate B^. at the missing point i.e., B , , /^^ /?~T77fr 
Along with the value oi B„ , , , /?,, we can now compute B, , ^, , , i r^, and in this manner we can iterate and 

6 2+1/V2' ^ (2+l/\/2) + (l/V2+l/\/2) 

obtain the entire solution. 

It should be noted that our interpolation based approach is geared towards finding only pre-classical solutions i.e. 
solutions that do not vary much between neighboring points. Restricting ourselves to pre-classical solutions in the 
context of this work is meaningful because our goal is to explore the semi-classical sector of LQC models. In fact, as 
pointed out earlier, lattice-refinement plays a significant role only in that regime. 

Complete solutions for the B^- relation for a few different values of the separation constant A appear in Figure [l] 
An immediate concern that one may have with our approach is the matter of the sensitivity of our results to the 
choice of an interpolation scheme i.e., the least-squares fit. In Figure |2] we plot the percent-difference, d between 
various solutions of the Br relation, obtained by varying the type of interpolation scheme. The results clearly indicate 
that our solution method is very robust (the percent-difference is at the 0.1% level). In the remaining portion of this 
article, we will continue to use a least-squares fit as our choice of the interpolation scheme. 

An alternate approach that is applicable to all ID lattice-refined relations [lOj lUj , suggests that we define a change 
of variable, f ~ It'^/^, one that can make the Bj- relation equi-spaced to linear order in S. This is because, in terms 



of this variable, B 



T±2SN^ 



Bf±2S and then the B^ relation \1\, to linear order in 5, becomes 



Br- 



-2(5 



Bf-2S~-{\+2)Br5/if. 



(3) 



This particular relation has been studied before at various places in the literature [5] and its solutions are well known. 
It appears in the separable Bianchi I LRS (LQC, without lattice refinement) model and the asymptotic form of its 
solution can be written analytically as: 



Br oc r 



-(A+2)/12 



OC T 



-(A+2)/8 



(4) 



In Figure [3] we plot the two solutions obtained using these different approaches together. The agreement is clearly 
seen to be excellent. 

We now move on towards a solution of the A^ relation ([2|. Note that our interpolation based method is immediately 
applicable to this relation as well. Indeed, this is the main merit of the technique; it is generically applicable to lattice- 
refined difference equations. Solutions for a few different values of the separation parameter A are plotted in Figure [4] 
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FIG. 3: We plot the solution of the B^ relation solved using the interpolation approach alongside the approximate analytic 
solution of the same relation. Here the solid-line is the analytic solution. There is an excellent agreement between the two 
approaches. 



Let us proceed towards finding an asymptotic, analytic solution by transforming the relation to an equi-spaced one, 
using the new variable: fi, = |/i'^/^. Using, A j^j^^^^-i w A^±4^s and also that fi, >> 6 the A^ relation can be simplified 
and written as an approximate ODE which takes the form shown below: 



(PA 
all 



' 444'' 



(A-2)/3 



dA 
dft, 



(5) 



This simplified equation, that can be assumed to be valid asymptotically can be solved easily yielding solutions in 
terms of modified Bessel functions: 



A^ (X m('+i)/4/(_,_i)/4(^^/2), m(^+i)/4X(,+i)/4(7/^/2) 



(6) 



It should be noted that these modified Bessel functions exhibit exponential behavior asymptotically, and it is related 
to this fact that the separable Schwarzschild model under consideration (with Nr oc y^frj and N^ ex V^H) exhibits 
generic instability |10j . In Figure p^ we plot this analytic solution alongside the numerically generated solution. The 
agreement is clearly seen to be excellent. 



B. Two-dimensional lattice-refined relations 



In the previous section, we solved the master Schwarzschild interior LQC relation for a simple intuitive choice 
for Nr and N^j^, one that makes the relation variable separable. We noted that the solutions to that version of the 
relation are unstable, which implies that the model cannot be viable as a physical model for the quantum geometry 
of the Schwarzschild black hole interior. This outcome is consistent with what was observed before [101 wherein, von 
Neumann numerical analysis of the same relation yielded the same unstable characteristics. Now, let us consider the 
alternate choice for N^ and iV^, the one that has improved stability properties: N^ ex \/\t\ and N^ ex /^/-\/|r[. It is 
very interesting that the notion of von Neumann stability that plays ordinarily a significant role in numerical analysis 
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FIG. 4: We plot the solution of the A^ relation for different values of the separation constant A, using the interpolation method 
as presented in the text. The dashed-line corresponds to a value of X = 10, the dotted line corresponds to X = 5 and the solid-line 
is for A = 1. 



of partial differential equations, provides an excellent mechanism for limiting the wide range of possible choices for 
Nr and Nf, [10]. 

With this choice of N^ and N^, the relation under consideration is not variable separable. This is because the 
step-size of any of the two independent variables depends on the other variable. In addition, to complicate matters 
further, it is not possible to find two independent quantities that would make it possible for the difference equation 
to be transformed into a uniformly spaced version [10', just like we did for the ID relations in the previous section. 
Thus, we have no choice but to solve this non-uniformly spaced, partial difference equation directly. Fortunately, 
the local interpolation method, that we have introduced in this article and detailed in the previous section is readily 
applicable to this challenging problem. 

Conceptually, our approach towards solving this 2D lattice-refined difference equation is identical to the one we 
took in ID. The relation connects various points on the lattice, which are shown clearly in Figure [6] In order to solve 
for the value '4'n+2SN'^ t+25n~^ ^'^ evolve the solution forward, we need to have values for the wave function that 
in general, have not been computed from previous iterations. However, these values are in close proximity to values 
that are known from before, therefore it is very natural to perform a local 2D interpolation using the values that are 
known, to compute the values that are needed. To demonstrate the success of this approach, we evolve a semiclassical 
wave-packet and show the results in Figure [7] We explain the evolution process in some detail below. 

We start by assuming a gaussian profile for the wave function ^, and evolve it using the interpolation approach 
that we developed in the previous section. A 2D grid for /z and r points is first set up by arbitrarily choosing some 
(/i, r) ordered pair with both /i and r being large positive integers. The grid implements the fact that /x advances 
as fj, + l/V^ while r advances as t -I- fx/^/r. Each point on the grid is labeled by two integers (m,n) which would 
represent the corresponding position of the point on an uniformly spaced grid. Note that although we don't have an 
uniformly spaced grid, these labels help with "book-keeping" i.e. identifying the points being used for the iterations. 
The first row of points that have a constant n label and are of the form (/x-t- l/\/^i ''") i-^- the set of points (relevant to 
the relation) for which r is a constant. All the other points evaluated on the grid are of the form (fj,+ l/^/T, r + fi/^/r). 
Each of these points act as the central anchor for the other points that enter the relation. For example, referring to 
Figure [6| it is the central point of the polygon that we evaluate while laying out the grid. All the points "connected" 
to it via the relation, can then be easily evaluated. 

Now that we have explained how the 2D grid is set up, let us turn the discussion towards evaluating and evolving 
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FIG. 5: We plot the logarithm of the solution of the A^ relation solved using the interpolation approach alongside the approximate 
analytic solution of the same relation. Here the solid-line is the analytic solution. There is an excellent agreement between the 
two approaches, especially for large ^. 



the wave function at these discrete points. As mentioned before, the initial wave function is assumed to be a gaussian 
wave-packet centered on the fj, axis. The values of the wave function at all the points marked as dark-circles in 
Figure [6] are evaluated using the gaussian wave-packet expression and this data is used to find the value of the wave 
function at the location (fi + 1/^/t, t + [ij \/t\ which is the white-circle in the figure. In this manner, we obtain the 
value of the function at points in the vicinity of (/i + 1/\/ti ''") for all /i. In other words (with respect to labels m, n), 
we compute the values of the wave function for all values of to, for a given n. Then we move on to the next n-step 
(recall that to and n are just integer labels for the points with respect to the starting point, so they advance by one 
along their respective axes). But, now we find ourselves in the same situation as we were in ID i.e. we do not have 
the data points we need to make the evolution progress further using the 2£) master relation. Thus, we utilize a 2D 
least-squares fit to compute these unknown values using the data points known from the previous n-step, that have 
the same value of to,. In other words, to find the value of the wave function that is needed by the ID relation, we find 
a local ID interpolation formula for points in the neighborhood of (m^n) using the data from the {m,n — 1) point 
and the points around it (see Figure p|. This yields an explicit scheme for progressing the evolution forward. 

The results from this approach are depicted in Figures Wj (for the Schwarzschild interior model) and Figure p^ (for 
the Bianchi I LRS model) . It is clear that the evolution of the wave-packet is very smooth and stable as one would 
expect in the semi-classical regime. It should be noted for both figures, we chose the domain to be such that /i > 4r, 
which happens to be the unstable region of the domain for the equi-spaced case [H]. Therefore, it is clearly evident 
that lattice-refinement has "cured" the problem of von Neumann instability in such anisotropic LQC models. This is 
completely consistent with the analytical stability analysis that was performed in |10j . The starting values we used for 
(/x, r) were (100000,25000) and we performed approximately 200 iterations in both directions (250 x 150 grid points). 
In the near future we plan on developing a high-performance numerical code based on this current proof- of- concept 
approach that will be able to compute on much larger grid sizes, in a very effective manner. Then we will be able 
to study the semi-classical wave-packet trajectories on a much larger domain and compare them with the recently 
published effective equations and effective dynamics [ T^ . 
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FIG. 6: The numerical "stencil" that is used for the evolution of a wave-packet in the lattice-refined LQC model of the 
Schwarzschild interior. This graphic depicts all the values of the wave function that are connected by the evolution differ- 
ence equation (Hamiltonian constraint). 



III. SUMMARY & DISCUSSION 



In this article, we have presented an effective numerical method for solving the Hamiltonian constraint in generic 
LQC models with lattice-refinement. The method is based on performing a /oca/ interpolation to obtain needed values 
from known values. 

We demonstrated the success of the method by solving the lattice-refined LQC model of the Schwarzschild interior 
geometry using two different refinement models. In the first refinement model, the 2D quantum evolution equation 
becomes variable-separable, reducing down to two ID lattice-refined equations that we then solved using our interpo- 
lation technique and also analytically, using asymptotic approximations. Using this model as a test case, we compared 
the accuracy and the robustness of our interpolation approach. Then we chose a refinement model that is inherently 
non-separable, and solved that in the context of the evolution of semi-classical wave-packets using an extension of the 
same interpolation approach in two dimensions. The results from our evolutions reinforced recently published results 
on the von Neumann stability analysis of the lattice-refined Schwarzschild relation [10]. 

One preliminary remark that we would like make is that, based on our current results we do not see any evidence 
of the kind of behavior depicted in |12j which was obtained from a study of the effective semi-classical dynamics. In 
particular, in [T^ the effective semi-classical Hamiltonian constraint is solved for the Schwarzschild interior model and 
those results suggest that: (a) In the equi-spaced model, there is a "bounce" and the quantum dynamics matches the 
solution to a different black hole (note that the "bounce" is not through the classical singularity, it actually avoids it!); 
and (b) In the lattice-refined case, one obtains an equilibrium solution asymptotically, thus avoiding the singularity 
(again). Results from the explicit evolution of wave-packets using the quantum Hamiltonian constraint appear to be 
qualitatively different. For example, in the equi-spaced case, past results [HI |9] suggest that the quantum evolution 
equation allows for an evolution right through the classical singularity, as opposed to avoiding it. In addition, our 
preliminary results from the lattice-refined case show no indication of an asymptotic equilibrium, again as suggested 
by |12j . However, we make these preliminary remarks with caution. There are many caveats associated to the results 
presented in |12j . as the authors themselves point out (basically, they include one quantum effect, but ignore others). 
In addition, the approach towards numerical evolutions of wave-packets in lattice-refined LQC that we are presenting 
in this article is just a proof-of-concept demonstration, and is currently lacking any serious analysis and study. We 
will attempt to explore this apparent "discrepancy" in the near future. 




FIG. 7: Stable evolution of a gaussian wave-packet in the lattice-refined LQC model of the Schwarzschild interior, based on the 
local 2D interpolation method as presented in this article. The scales of the axis depicted above are based on (m, n) and are 
arbitrary. The accompanying text has more detail. 
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FIG. 8: Stable evolution of a gaussian wave-packet in the lattice-refined LQC model of Bianchi I LRS, based on the local 2D 
interpolation method as presented in this article. The scales on the axis depicted above are based on (m, n) and are arbitrary. 
The accompanying text has more detail. 



